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ABSTRACT 

We  present  a  systematic  method  for  constructing  boundary  conditions  (numerical  and  physical) 
of  the  required  accuracy,  for  compact  (Pade-like)  high-order  finite-difference  schemes  for  hyperbolic 
systems.  First  a  proper  summation-by-parts  formula  is  found  for  the  approximate  derivative.  A 
“simultaneous  approximation  term”  (SAT)  is  then  introduced  to  treat  the  boundary  conditions. 
This  procedure  leads  to  time-stable  schemes  even  in  the  system  case.  An  explicit  construction  of 
the  fourth-order  compact  case  is  given.  Numerical  studies  are  presented  to  verify  the  efficacy  of  the 
approach. 
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Introduction 


Emphasis  on  the  long-time  numerical  integration  of  the  fluid  mechanics  equations  lias  increased 
in  recent  years.  As  a  result,  high-order  spatially  accurate  schemes  are  favored,  because  of  their  lower 
phase  error.  Such  schemes,  although  they  are  stable  in  the  classical  sense  (Lax  and  G-K-S  stability), 
may  exhibit  a  non-physical  growth  in  time.  For  a  fixed  time  T,  these  schemes  converge  as  the  mesh 
size  Air  — >  0.  However,  from  a  practical  point  of  view,  in  order  to  achieve  reasonable  accuracy  for 
large  T,  meshes  much  too  fine  for  the  computers  available  in  the  foreseeable  future  are  required. 
Since  long-time  integrations  are  encountered  in  present  day  computations,  it  is  important  to  devise 
schemes  which  are  not  only  classically  stable  but  also  time-stable.  Specifically,  they  do  not  allow  a 
growth  in  time  that  is  not  called  for  by  the  differential  equations. 

To  retain  the  formal  accuracy  of  a  high-order  scheme,  boundary  closures  must  be  accomplished 
with  accuracies  that  are  at  most  one  order  less  than  the  interior  scheme  [1].  For  the  scalar  explicit 
central-differencing  case,  Kreiss  and  Scherer  [2]  have  presented  a  method  for  constructing  a  boundary 
condition  of  accuracy  one  order  less  than  the  inner  scheme  such  that  a  generalized  summation-by-parts 
property  of  the  differential  equation  is  preserved.  Strand  [3]  has  used  their  approach  1o  construct 
in  the  scalar  case,  fourth-  and  sixth-order  central-differencing  schemes  with  boundary  closures  of 
the  appropriate  order  such  that  the  resulting  expression  for  '  V  derivative  satisfies  the  summation- 
by-parts  property.  Recent  attempts  to  utilize  these  boundary  closures  to  numerically  solve  a  2  x  2 
hyperbolic  system  have  shown  that,  in  certain  cases,  an  unwarranted  growth  in  time  still  results. 

In  reference  [4],  the  stability  characteristic  of  various  compact  fourth-  and  sixth-order  spatial 
operators  were  assessed  using  the  theory  of  Gustafsson,  Kreiss  and  Sundstrom  (G-K-S)  [5]  for  the 
semidiscrete  initial-boundary-value-problem  (IBVP).  This  study  showed  that  many  of  the  higher 
order  schemes  that  are  G-K-S  stable  are  not  time  stable.  It  was  concluded  that,  in  practical  calcula¬ 
tions,  only  those  schemes  which  satisfied  both  definitions  of  stability  were  of  any  usefulness  for  long 
time  integrations.  Of  practical  importance  was  a  new  sixth-order  scheme  with  fifth-order  boundary 
conditions  which  was  shown  to  be  G-K-S  and  time-stable.  Recently,  however,  it  has  been  found  that 
most  of  the  high-order  schemes  that  were  time-stable  in  the  scalar  case,  exhibited  time  divergence 
when  applied  to  a  2  x  2  system. 

In  this  paper,  we  outline  a  systematic  procedure  for  designing  time-stable,  as  well  as  G-K-S 
stable  schemes  of  high-order  accuracy.  The  new  schemes  are  guaranteed  to  be  t  ime-stable  for  any 
hyperbolic  system  (as  long  as  the  system  has  a  bounded  energy).  The  first  step  in  this  procedure  is 
to  construct  an  approximation  to  the  first  derivative  (internal  plus  boundary  points)  that  admits  a 
summation-by-parts  formula.  We  l  dy  on  the  work  of  Strand  [3]  for  high-order  explicit,  formulations. 
For  high-order  compact  schemes,  we  derive  a  new  methodology  for  construction  of  such  schemes. 


Appendix  I  includes  an  exposition  of  the  methodology,  and  a  detailed  example  of  the  fourth-order 
compact  central  difference  scheme  with  third-order  boundary  closures.  In  section  1 ,  we  discuss  a  scalar 
hyperbolic  equation.  We  show  that  in  general  a  summation-by-parts  formula  does  not  guarantee  time 
stability.  However,  we  introduce  a  new  procedure  for  imposing  boundary  conditions  (simultaneous 
approximation  term,  (SAT)),  that  solves  a  linear  combination  of  the  boundary  conditions  and  the 
differential  equations  near  the  boundary.  This  technique  is  an  extension  of  the  techniques  used  in 
reference  [6]  to  stabilize  the  pseudo-spectral  Chebychev  collocation  method.  It  is  shown  that  if  the 
approximation  of  the  derivative  operator  admits  a  summation-by-parts  formula  then  the  SAT  method 
is  stable  in  the  classical  sense  and  is  also  time-stable. 

In  section  2  we  discuss  the  implementation  of  the  SAT  method  to  systems  of  hyperbolic  equations. 
We  show  that  also  in  the  system  case,  time  stability  (as  well  as  Lax  stability)  is  assured  by  having  a 
summation-by-parts  property  for  the  numerical  derivative  operator,  provided  that  the  SAT  method 
is  utilized. 

In  section  3  we  present  numerical  results  that  confirm  the  efficacy  of  the  SAT  procedure  even  in 
the  cases  where  previous  attempts  could  not  attain  time  stability.  It  is  shown  that  the  theoretical 
predictions  for  the  time  stability  of  the  SAT  method  are  realized  in  practice  for  both  the  scalar 
hyperbolic  case  and  the  2x2  hyperbolic  system.  Finally,  an  optimization  of  the  parameter  r  (which 
arises  in  the  SAT  procedure)  is  performed,  with  regard  to  efficiency  and  accuracy. 

1.  The  Scalar  Case 


We  consider  the  scalar  hyperbolic  equation 

du_  du 
dt  dx 


0  <  *  <  1 


for  which  there  exists  the  energy  rate 

—  /  u2(x,  t)dx  =  A(u2(l,  t)  —  u2(0,  t)) 
dt  Jo 

For  positive  A,  we  have  the  boundary  condition 


(1) 


«(M)  =  g(t) 

We  denote  by  u  a  vector  of  the  unknowns  (u0(t),U](t),  ...u^(t))  which  corresponds  to  grid  points 
*o(=  0),xj,  ...xN(=  1). 


In  this  work,  we  deal  primarily  with  compact  schemes  for  the  discretization  of  the  spatial  operator 
Fcr  a  compact  spatial  operator,  the  approximation  to  the  first  derivative  can  be  written  as 


du 

P  — =  QU 
dx 


(2) 


where  P  and  Q  are  ( N  +  1)  x  ( N  +  1)  matrices.  We  further  assume  that: 


Assumption  I 

(i)  Equation  (2)  is  accurate  to  order  m.  Specifically,  if  we  denote  by  v  the  vector  (u(xo,  t ), v(x^, 
where  v(x,t)  G  Cm  and  Xj  =  j Ax  =  -fa  ,  and  by  vx  the  values  of  ((|^)0,  •••,  (|^)Af)r  then 

Pvx  -Qv=  PTe 

where  the  truncation  error  Te  satisfies 

|Te|  =  0(  Ax)m 

(ii)  The  matrix  P  has  a  simple  structure  (preferably  tridiagonal)  and  is  easily  invertible. 

(iii)  There  exists  a  matrix  H,  and  positive  constants  pi,  y,2  independent  of  N  such  that 

fill  <  HP  <  fi^I 

specifically,  HP  is  a  symmetric  positive  definite  matrix. 

(iv)  There  exists  a  matrix  G  —  H  Q  such  that  G  +  GT  has  only  two  elements:  g0fi  and  giv^. 

In  general  we  require  g0,o  <  0  <  #/v,yv- 

Assumptions  1  and  2  are  common  to  any  useful  compact  scheme.  Assumptions  3  and  4  are  specific 
to  the  summation-by-parts  requirement  for  the  spatial  operator. 

Equation  (1)  is  now  semi-discretized  using  formula  (2)  to  yield 

~  =  AP-'Qu  (3) 

Note  that  assumptions  3  and  4  from  above  admit  a  summation-by-parts  formula  in  the  sense  that 

dE  2  2  f*\ 

=  9o,ou0  +  gN,NUN  (4) 

where 

£(()  =  i(u(<),tf/>u«)  (5) 

In  Appendix  I  we  show  how  to  construct  a  fourth-order  compact  scheme  that  satisfy  Assumption 
1  and  therefore  (4). 

Interestingly,  equations  (4)  and  (5)  were  obtained  without  imposing  the  boundary  conditions.  We 
will  use  the  summation-by-parts  property  defined  in  equations  (4)  and  (5)  to  construct  a  scheme 
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that  admits  a  decreasing  energy  norm  when  the  boundary  condition  is  imposed.  Note  that  the 
way  in  which  the  boundary  condition  is  imposed  is  important  for  numerical  stability.  The  most 
common  procedure  of  imposing  the  boundary  conditions  (A  >  0  ).  is  to  use  equation  (3)  to  update 
the  unknowns  u0,  ..u/v,  followed  by  overwriting  u y  =  y{t).  This  procedure  accounts  for  the  fact  that 
in  a  general  hyperbolic  system  the  precise  location  for  each  boundary  condition  is  not  known  until 
after  a  characteristic  decomposition  is  performed  at  all  boundaries.  This  procedure  (particularly  if 
//  is  a  nontrivial  matrix),  may  not  yield  the  estimate  (4)  with  u \  replaced  by  </(/).  In  short,  the 
imposition  of  certain  boundary  treatments  may  ruin  the  structure  of  the  summation  norm,  which 
results  in  a  numerical  scheme  that  is  not  time-stable. 

A  simple  counter-example  is  presented  which  demonstrates  the  necessity  of  careful  boundary 
implementation.  Consider  the  scalar  equation  ut  =  ux  with  the  boundary  condition  a \  =  </(t).  The 
semi-discretization  in  the  absence  of  boundary  conditions  becomes  ut  =  A  it.  where  .4  =  P-1  Q.  As 
described  earlier,  once  the  matrix  A  is  formed,  the  boundary  conditions  are  imposed.  This  has  the 
effect  of  pre-multiplying  the  matrix  A  by  the  boundary  matrix  D.  Without  loss  of  generality,  we  use 
the  boundary  condition  g(t )  =  0  in  this  problem;  the  resulting  boundary  operator  is  the  matrix 


D: , 


1  0  0 
0  1  0 
0  0  0 


For  time  stability,  the  resulting  matrix  ,4t  =  D  P  1  Q,  rather  than  the  matrix  4  must  exhibit  a 
summation- by-part  norm. 

for  simplicity,  we  discretize  the  domain  into  two  even  intervals,  such  that  the  discrete  solution 
r 

vector  is  (u0,  U\,  112)  .  T  he  boundary  condition  is  imposed  at  u2.  A  first-order  discretization  that 
satisfies  the  summation-by-parts  energy  norm  is 


77/48 

( — 19)/ 12 

( — 43)/ 48  ' 

'  (—25)/ 10  4  (—39)/ 10  ■ 

p»  = 

(  —  19)/ 12 

32/3 

(  —  13)/ 12 

;  <?:.  - 

-4  0  4 

( — 43)/  48 

(  I 3)/ 1 2 

53/48 

39/16  -4  25/16 

Note  that  the  matrices  P  and  Q  satisfy  P:i  =  Pj  and  Q:i  =  except  for  7 ,M,  and  1/2,2.  In 

this  example,  the  matrix  //  is  the  identity  matrix.  The  characteristic  equation  for  the  P:i  matrix  is 
—  19‘2A3 +  2568A2  —  5026 A  +  501  =  0.  The  symmetry  of  Pt  and  the  alternating  signs  of  the  respective 
terms  in  the  characteristic  polynomial  guarantee  the  positive  definiteness  of  P(.  The  discretization 
operator  4,j  =  Pp 1  Q:\  can  be  written  as 


1 


11/1002  ( — 512)/ 50 1  1013/1002  * 

A3  =  (— 55)/334  ( —  1 12)/ 167  279/334 

.  2059/1002  (—2560)/. 501  3061/1002 

All  the  requirements  of  the  summatiou-by-parts  energy  norm  are  satisfied  by  this  discretization,  and 
a  precise  energy  norm  exists  in  the  absence  of  boundary  conditions. 

The  combined  operator  A\  =  D3  A3  becomes 

11/1002  ( — 512) /501  1013/1002 

A\  =  (— 55)/334  (-112)/ 167  279/334 

0  0  0 

for  which  the  characteristic  polynomial  is  —  1002A3  —  66 1  A2  + 1 76A  =  0.  The  roots  of  the  characteristic 
polynomial  are  A  =  —0.86317...,  and  A  =  0.20349...,  respectively.  The  numerical  solution  will  grow 
in  time  as  a  result  of  the  eigenvalue  in  the  right  half  of  the  complex  plane  (RH-P)  and  will  not  be 
time-stable. 

As  demonstrated  by  the  previous  counter-example,  a  spatial  operator  which  satisfies  the  summation- 
by-parts  energy  norm  may  not  be  time-stable.  Many  of  the  high-order  schemes  that  satisfy  the  sum¬ 
mation  property  are  time-stable  for  the  scalar  case.  A  notable  exception  is  the  sixth-order  explicit 
scheme  with  fifth-order  boundary  conditions  reported  in  the  work  of  Strand  [3].  (See  Appendix  II  for 
details  of  this  scheme.)  For  this  sixth-order  scheme,  time  stability  can  be  guaranteed  only  if  the  last 
row  and  column  of  the  matrices  HP  and  HQ  are  removed  before  matrix  inversion  and  multiplication 
are  performed. 

The  underlying  reason  for  the  growth  in  time  is  the  imposition  of  the  boundary  condition  operator, 
which  has  an  effect  on  the  structure  of  the  norm  matrix  P  in  ut  —  D  P_1  Q.  Specifically,  D  P~l 
destroys  the  structure  of  the  norm  P.  In  the  scalar  case,  this  problem  can  be  eliminated  in  certain 
circumstances.  For  instance,  if  the  matrix  P  is  a  restricted  full  norm,  then  D  P-'  still  produces  a 
useful  norm  by  eliminating  the  zero  element.  A  restricted  full  norm  is  defined  where  the  diagonal 
is  the  only  nonzero  element  in  the  first  (or  last)  row  and  column  of  the  matrix  P  (See  Strand  [3]). 

A  special  case  of  the  restricted  full  norm  is  the  diagonal  case,  which  is  of  some  practical  interest. 
Unfortunately,  even  for  cases  where  P  is  a  restricted  full  norm,  stability  cannot  be  generalized  to  the 
case  of  a  hyperbolic  system.  An  alternative  means  of  imposing  boundary  conditions  must  be  found 
for  these  cases. 

At  this  point,  we  introduce  the  SAT  methodology  for  boundary  implementation.  We  show  in 
the  following  text  that  the  SAT  method  leads  not  only  to  stability  but  also  to  time  stability  for  the 
scalar  wave  equations,  and  this  property  applies  to  arbitrary  hyperbolic  systems.  The  SAT  method 


involves  the  indirect  imposition  of  the  physical  boundary  conditions.  This  is  accomplished  by  adding 
a  term  to  the  derivative  operator,  which  is  proportional  to  the  difference  between  the  discrete  value 
un  and  the  boundary  term  g(t).  Thus,  we  propose  the  discretization, 

=  \Qu-TXgNtNS(uN  -  g{t))  (6) 

where 

S  =  0, ...,  0, 1)T  (7) 

Contrary  to  the  common  practice  of  satisfying  the  boundary  condition  directly  by  imposing  u h  =  g(t), 
the  SAT  method  involves  solving  a  derivative  equation  everywhere ,  including  the  boundary  points. 
The  extra  term  which  is  added  accounts  for  the  boundary  information  to  within  the  accuracy  of  the 
original  discretization.  Note  that  the  SAT  is  added  not  only  to  the  boundary  equation  but  to  other 
points  depending  on  the  structure  of  the  vector  S,  (which  is  the  last  column  of  the  matrix  H~‘). 
The  extra  SAT  term  does  not  alter  the  accuracy  of  the  scheme,  since  the  SAT  term  vanishes  upon 
substitution  of  the  analytic  solution. 


We  now  demonstrate  that  the  SAT  method  yields  a  Lax  stable  and  time-stable  scheme.  For  the 
time  stability  analysis,  we  take  g(t)  =  0.  We  pre-multiply  equation  (6)  by  H  and  use  equation  (7) 
to  obtain 


HP^  =  XHQu  -  tA<7jv,/v(0,  0,  ...,0,  l)rUAf 
at 

We  now  define  the  energy  E(t)  as  in  equation  (5)  to  get 

d  E  ( t )  2  2  2 

— ^ —  =  9o,ou0  +  9n,nuiw  ~  t9n,nun 

With  </o,o  <  0  <  <7amv,  we  can  immediately  state  the  following  theorem. 


(8) 

(9) 


Theorem  1.1: 


The  SAT  method  presented  in  equation  (6)  is  both  stable  and  time-stable  if 

r  >  1 


(10) 


In  addition  to  proving  the  stability  of  the  SAT  scheme  defined  in  equation  (6),  we  must  show  that 
the  procedure  preserves  the  order  of  accuracy  m  of  the  spatial  operator.  This  is  accomplished  by  a 
direct  convergence  proof  showing  that  the  SAT  term  indeed  preserves  the  spatial  order  of  accuracy. 

Denote  by  v  the  vector  (u(xo,t),  ...,u(xN,t))T  ,  i-e.  the  values  of  the  true  solution  of  (1)  at  the 
grid  points.  Combining  the  accuracy  condition  found  in  Assumption  /  with  equation  (6)  we  have 

d\ 

F—  =  \Qv  -  r\gN<NS  [ufarjv,  <)  ~  g{t)]  +  P Te 
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(11) 


Note  that  u(xn,  t)  —  g(t )  =  u(l,  t)  —  g(t)  =  0.  Now  define 


Cj{t)  =  u(xj,t )  -  u3{t) 


where  Uj(t)  solves  (6)  ,  to  obtain 


P7t=  ~  T^9N’NSeN  +  -PTe 


(12) 


where  Te  is  the  truncation  error  defined  in  Assumption  1 .  We  now  use  the  energy  estimate  presented 
in  (9)  to  obtain 

d{e,  HPe ) 


dt 


<  (e,  HPTe) 


and  the  inequality 


to  obtain 


(e,  HPTt)  <  y/(c,  TTpT)J{T^HPT)) 


<  yjiT.'HPT,) 


(13) 


By  assumption  I  ,the  truncation  error  is  of  order  m,  and  we  get 

J{t~HPt)  <  0( Ax)m 

which  proves  the  convergence  of  the  scheme. 


In  conclusion,  a  precise  means  is  now  available  for  the  scalar  case  to  impose  boundary  condi¬ 
tions  that  are  guaranteed  to  be  time  stable,  and  that  preserve  the  formal  accuracy  of  the  original 
discretization. 


2.  The  Hyperbolic  System 


In  this  section,  we  explain  how  to  use  the  SAT  method  for  systems  of  hyperbolic  equations  and 
show  that  the  resulting  scheme  satisfies  an  energy  estimate  similar  to  the  one  obtained  for  the  scalar 
differential  equation.  First  the  system  of  differential  equations  is  described. 

Let  u1  and  u11  be  the  two  function-valued  vectors 

u1  =  (lt(l)(x,f u(k)(x,t)) 
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(14) 


U{T\x,t)) 


uH  =  (/+’),.., 

that  solve  the  system  of  differential  equations 

du1  ,  du1 

- =  A  - 

dt  dx 


du11  A  u  du11 

— - —  =  A  - 

dt  dx 

where  A7  and  A11  are  diagonal  matrices  of  the  form 

A7  =  diag(\i,...,\k) 

A"  =  diag{\k+u-~,Ar) 

In  order  to  impose  the  boundary  conditions  we  assume  that 

A,  >  A2  >  ...  >  A*  >  0  >  Afc+1  >  ...  >  Ar 

For  this  case,  a  well-posed  set  of  boundary  conditions  is  given  by 

ur(M)  =  R\ill{l,t)  -f  gX(t) 


(15) 


(16) 


un(0,0  =  Lu^OJ)  +  gn(t) 

Where 

gI(t)  =  (g(1»(t),...,g(k>(t)) 

and, 

g!I(t)  =  (g"t+,>(t),....g(r»(t)) 


(IT) 


In  equation  (17),  the  matrix  R  has  k  rows  and  r  —  k  columns,  while  the  matrix  L  has  r  —  k  rows 
and  k  columns.  Without  loss  of  generality,  for  the  stability  analysis  we  will  assume  that  both  g^t) 
and  gXI(t)  vanish. 

Equation  (17)  is  well-posed  for  any  L  and  R.  However,  to  guarantee  no  growth  in  time  some 
conditions  must  be  imposed  on  the  matrices  L  and  R.  These  conditions  are 
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Condition  I: 


\L\\R\<\  (18) 

where  the  matrix  norm  is  defined  by 

\A\  =p{ATA)i 

and  p(A)  is  the  spectral  radius  of  A. 

The  Continuous  case 


It  is  instructive  to  establish  and  prove  an  energy  estimate  for  the  continuous  hyperbolic  system 
although  such  a  proof  is  well  known.  The  same  basic  steps  that  are  used  in  the  continuous  proof  will 
be  used  later  in  the  text  to  prove  the  energy  estimate  resulting  from  the  semi-discrete  hyperbolic 
system. 

Condition  I  is  a  sufficient  condition  for  the  solution  of  equation  (15)  to  be  bounded  in  time.  In 
fact  one  can  state 


Theorem  2.1: 


Let  and  u11^,*)  be  the  solution  of  equation  (15)  with  the  boundary  conditions  (17). 

Recall  that  we  take  g1  =  g11  =  0.  Suppose  that  L  and  R  in  equation  (17)  satisfy  Condition  I.  Define 
an  inner  product 


and  an  energy  function  E(t) 


(w,v)  =  [ 
Jo 


w(x,  t)v(x,  t)dx 


(19) 


£(<)  =  £  ^’V0)  +  t  }£[<““’ 

«'=i  A«  »=fe+i 


,(*>1 


then  the  time  rate  of  the  energy  function  satisfies 


dE 

dt 


<  0 


(20) 


(21) 


Proof 


We  start  by  differentiating  equation  (19)  with  respect  to  t  to  obtain 

d(u^\  u(0) 


dt 


=  2  f  «0)u|0 

Jo 


dx 
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Using  equation  (15)  we  obtain 


so  that 


dt 


=  2  /’  u^XiU^dx 
Jo 


d(u^Ku^)  /-v  2  ,  1  2 

- — - =  A,(u(,)(l,t)  -u(,)(0,f))  (22) 

Differentiating  equation  (20)  and  substituting  equation  (22)  we  obtain  the  energy  rate  for  the  system 
as 

—  =  £|i|(u'i|(l,()2-U|i,(0,i)2)-  £  |fl|(ul'l(l,i)2-uW(0,i)2)  (23) 

«=1  t=fc+l 

relating  the  time  rate  of  change  of  the  energy  function  to  the  energy  that  crosses  the  boundaries. 
Note  the  change  of  sign  in  the  second  term  which  results  from  the  negative  sign  of  the  eigenvalues 
X i  for  k  <  i.  We  must  now  quantify  the  magnitude  of  the  boundary  terms  in  equation  (23). 

Replacing  the  sums  in  equation  (23)  with  the  vector  operations 

£*(i>(M)J  =  uI(l,<)TuI(l,l) 


£  u(,*(0,  t)2  =  uII(0,  i)Tun(0,  () 

«=fc+l 

we  can  now  make  use  of  the  boundary  conditions  in  equation  (17)  to  obtain 

uI(l,<)TuI(l,f)  =  un(  1 ,  t)T  RT  Run(l ,  t) 


un(0,t)Tun(0,0  =  u!(0,  f)TLTLuI(0,  t) 

Substituting  the  equations  (24)  and  (25)  into  (23)  we  obtain 

?  =  uII(l.()T{flTfi|i|  -  |fl|)un(l,i)  +  {LT L\R\  -  |i|}u"(0,() 


Because  Condition  I  ensures  that 


RtR\L\  -  \R\I  <  0 


LtL\R\  -  I\L\  <  0 


and 


equation  (21)  is  established.  Therefore  the  continuous  energy  function  E{t)  is  bounded  in  time.  This 
completes  the  proof  of  Theorem  (2.1). 

The  Semi-discrete  Case 


We  are  ready  now  to  discuss  the  implementation  of  the  SAT  technique  for  the  system  in  equation 
(15)  with  the  boundary  condition  given  in  equation  (17).  As  in  section  1  we  denote  by  u‘  a  vector  of 
unknowns  (uq\ Uj\  which  correspond  to  the  grid  points  xo(=  0), x\,  ...x^{=  1).  We  assume 

that  we  have  matrices  P  ,  Q  and  H  such  that  the  scalar  case  admits  a  summation-by-parts  energy 
norm  given  in  section  1.  The  SAT  discretization  of  equations  (15)  -  (17)  is  chosen  as 

=  A ,gu‘  -  gN,N\tTS^(u^  -  (flu11)*?  -  sW)  1  <  i  <  k 


(27) 


P^jT  =  kQu'  -  5,o.oA,r6’(,)(4')  -  (£ur)o ]  -  9('])  k+l <i <r 

where  r  is  a  stabilizing  factor  to  be  determined  later.  As  in  the  scalar  case,  we  choose  S ^  to  be  one 
of  the  vectors 


=  H~\ 0,0,...,0,1)T  1  <i<k 


(28) 


S{t)  =  H~l(l,  0,...,0,0)T  *+l<i<r 

We  recall  from  the  scalar  case  that  HP  is  symmetric  positive  definite  and  HQ  is  skew  symmetric 
except  for  the  terms  </o,o  =  {HQ) o,o  <  0  and  <7n,n  =  {HQ)w}n  >  0.  Thus  equation  (27)  is  well 
defined. 


Before  proving  the  stability  (and  time  stability)  of  the  SAT  method  in  equation  (27),  we  would  like 
to  comment  on  the  role  of  the  matrix  H .  Explicit  knowledge  of  H  is  required  for  the  implementation 
of  the  SAT  method,  specifically  the  knowledge  of  <70,0  and  as  well  as  the  vectors  5'*)  are  needed 
to  implement  equation  (27).  Thus  H  is  not  only  a  theoretical  tool  (as  in  reference  [2])  but  is  also  of 
practical  importance. 

We  are  now  ready  for  the  stability  proof  of  the  SAT  method  in  equation  (27).  The  proof  is 
analogous  to  Theorem  2.1  with  the  continuous  integrals  replaced  by  discrete  sums.  The  scalar 
product  is  defined,  analogous  to  equation  (19),  as 

(u‘,u*)  =  X>!')u/(,) 

1=0 
11 


(29) 


A  different  scalar  product  to  be  used  later,  analogous  to  equations  (24),  is 


u1,  u1]™  =  £  U 


(•V0 

in  m 


[un,unl„,  =  E 

i=k+ 1 

for  m  =  0,  N. 

Theorem  (2.2) 

Let  the  SAT  method  defined  by  equation  (27)  satisfy  Assumption  1,  for  the  discretization  of  the 
hyperbolic  system  defined  in  equation  (15)  with  boundary  conditions  (17),  (w  ith  g*(t)  =  gn(t)  =  0). 
Then  the  discretization  is  both  stable  and  time-stable  provided  that 

2  -  2^/l-|fl||i|  2  +  2jl-|fi||£| 

l«l|£|  -T“  l*l|£|  (31) 


Moreover,  let  the  discrete  energy  be  defined  as 

*wo  =  Eyl(ui,//Pu')+  E 

*'=i  A'  i=k+ 1  lA,l 

where  the  scalar  product  (u\u‘)  is  defined  in  equation  (29).  Then 

dEN(t) 


Proof 


As  in  theorem  2.1  we  differentiate  the  scalar  product  (iT,  H Pu')  and  use  equation  (27)  to  obtain 
4(u\  HPu')  =  A,(u',  HQu*)  -  9n,nKt(u$  -  (Ruu)$)(u\  HS <‘>)  1  <i<k 


~{u\  H  Pu')  =  A,(u\  HQu')  -  go,oXT(u(0i]  -  (Rul)^])(ul,  HS{'])  k+l<i<r 


We  now  use  the  definition  of  from  equation  (28)  and  the  properties  of  HQ  from  Assumption 
I  to  obtain 

4(u',W/=u')  = 
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ito.oAitf)2  +  K9n,n(un])2  ~  K9n,nt(un])2  +  Xt9N,NTuy(RuU)^\  1  <  i  <  k 


(34) 


-<7o,o|A,|(uq,))2  -  \Xi\gN,N(u{N  )2  +  |At|^o,o7-(4!J)2  -  \Xl\go,oTii[o,(Lu1)(o\  k  +  1  <  i  <  r 

Note  that  in  equation  (34)  we  used  the  fact  that  the  A,  are  negative  for  A. •  +  1  <  i  <  r.  We  must 
now  quantify  the  magnitude  of  the  boundary  terms  in  equations  (34).  If  the  sums  in  equations  (34) 
are  replaced  with  the  vector  operations  defined  in  equations  (30)  we  get  an  estimate  for  the  discrete 
energy  rate 


,(d\2 


-(u-,HPu.)  = 
(<)\2 


=  l^l^o.ofu1^1]©  +  |I|fl-/v,/v(l  -  r  )[uI,  u1]^  +  (Ll^v./vfu1,  Ru11 


+  |/?|(r  -  l)flf0,o[uII,uII]o  -  |H|3^,Ar[un,un]N  -  g0fi\R\T[uu,  LuJ}0 
Substituting  the  estimates 

[u\Run]N  <  |uIljV|/2||uII|iv 
[u11,  Lux]0  <  |uII|o|I»||uI|0 

where 

l^lm  =  y^U^U1]™. 

into  equation  (35),  and  collecting  like  terms  yields 
dE^(t) 


N 


<  ~9n,n{\L\{t  -  l)^1^  -  r|L||rt||u1|iV|uu|jv  +  |/?||uu|^} 


Ii  u,IIi 


,H|2 


dt 


+£o,o{|fl|(r  -  l)|un|2  -  7'|i!i^||u1|o|u11|o  +  |Z/)|uA|g} 


Ii  1..II1 


,I|2l 


For  to  be  negative  we  require  each  curly  bracket  to  be  positive.  Thus  we  need 

\L\(t  -  l)^1^  -  r|L||.R||uI|/V|uII|;v  +  |*||unft  >  0 

and  also 

\R\iT  -  i)!uIIlo  -  T!Lil^iluI|o|un!o  +  l^llu1^  >  0 

Both  inequalities  are  satisfied  if 

|fi||i|rJ<4(r-l) 

and  this  is  equivalent  to  equation  (31).  Thus,  the  proof  is  established. 
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(35) 


(36) 


3.  Results 


Conventional  Boundary  Conditions 


Three  high-order  spatial  discretizations  (two  explicit  and  one  compact)  are  the  focus  of  the  results 
section:  the  fourth-order  explicit  scheme  with  third-order  boundary  conditions,  the  fourth-order 
compact  scheme  with  third-order  boundary  conditions,  and  the  sixth-order  explicit  scheme  with 
fifth-order  boundary  conditions.  All  satisfy  the  summation-by-parts  requirement  in  the  absence  of 
physical  boundary  conditions.  The  fourth-order  explicit  scheme  is  reported  elsewhere  (see  [3]  or  [7] 
for  specific  details)  and  will  not  be  derived  here.  The  fourth-order  compact  scheme  is  new,  and 
a  systematic  procedure  for  deriving  both  it  and  other  compact  high-order  schemes  is  presented  in 
Appendix  I.  The  sixth-order  explicit  scheme  was  first  reported  in  reference  [3],  but  is  also  included 
in  Appendix  II. 

First  we  demonstrate  that  all  three  schemes  behave  in  accordance  with  their  respective  order 
properties.  We  then  comment  with  regard  to  the  sixth-order  explicit  scheme,  that  satisfying  the 
summation-by-parts  energy  norm  is  not  sufficient  for  time  stability. 

The  model  problem  used  to  test  the  three  schemes  is  the  scalar  hyperbolic  equation 

du  du 

S+S  =  °,°<I<U>»  (37) 


u(0, t)  =  sin27r(— f),  t>0 


u(a:,0)  =  sin27r(ar),  0  <  x  <  1,  (39) 

The  exact  solution  is 

u(x,  t)  =  sin  27t(x  —  f),  0  <  x  <  1,  t>  0  (40) 

For  all  calculations,  the  time  discretization  used  was  a  fourth-order  Runge-Kutta  (R-K)  method 
with  the  time  step  small  enough  such  that  the  temporal  errors  are  much  smaller  than  the  spatial 
truncation  error.  In  all  cases,  the  boundary  condition  was  implemented  at  the  end  of  each  R-K  stage 
by  overwriting  the  value  of  the  solution  at  the  boundary  point. 

Table  I  shows  a  grid  refinement  study  performed  on  equation  (37)  for  all  three  spatial  discretiza¬ 
tions.  Both  the  absolute  (log  L2)  error  at  a  fixed  time  T  and  the  convergence  rate  between  two 
successive  grid  densities  are  plotted. 
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(fourth 

explicit) 

(fourth 

compact) 

(sixth 

explicit) 

Grid 

log  L2 

Rate 

Rate 

log  L2 

Rate 

Bi 

-0.501 

-1.418 

EH 

-2.080 

8.96 

-2.133 

4.06 

l 

1.88 

41 

-2.607 

4.22 

-2.627 

3.95 

7.29 

61 

-3.329 

4.10 

-3.316 

3.91 

-1.302 

8.17 

81 

-3.832 

4.03 

-3.806 

3.92 

-1.798 

3.96 

Table  I:  Grid  convergence  of  three  high-order  schemes  on  ut  +  ux  =  0. 


This  refinement  study  suggests  that  all  three  schemes  are  Lax  stable  (the  exact  solution  is  approached 
at  a  fixed  time  T  as  mesh  is  refined)  and  grid  converge  consistent  with  each  respective  theoretical 
rate.  The  convergence  rates  for  both  of  the  fourth-order  schemes  asymptote  to  the  theoretical  value 
of  4.  The  convergence  rate  of  the  sixth-order  explicit  scheme  is  sporadic  but  is  approximately  6 
(5.28  for  the  interval  between  21  and  81  points).  This  spurious  behavior  results  from  the  exponential 
divergence  of  the  solution  for  long  times  T .  At  T  =  70,  the  absolute  error  of  the  two  fourth-order 
schemes  is  comparable;  however,  that  of  the  sixth-order  scheme  is  two  to  three  orders  of  magnitude 
larger. 

These  numerical  results  indicate  that  the  two  fourth-order  schemes  are  time-stable;  the  sixth- 
order  scheme  is  not.  Nothing  in  the  definition  of  Lax  stability  precludes  exponential  divergence  of 
the  solution  for  long  times  T  as  long  as  the  divergence  rate  is  bounded  independently  of  the  grid 
used.  (See  reference  [4].)  The  numerical  divergence  of  the  solution  results  from  a  spatial  operator 
matrix  which  has  an  eigenvalue  with  a  positive  real  part  (an  RH-P  eigenvalue).  For  long  times  T, 
the  solution  is  dominated  by  this  eigenvalue. 

To  quantify  this  assertion,  a  comparison  is  presented  between  the  numerically  observed  divergence 
rate,  and  a  theoretical  prediction  from  eigenvalue  analysis.  By  assuming  that  the  numerical  error  can 
be  represented  as  ca r(t)  =  ejv(0)eaNt,  a  growth  rate  is  determined.  Similarly,  an  effective  growth 
rate  as  defined  by  easM At  =  \Gmax(At)\M ,  is  calculated  from  an  eigenvalue  determination.  (See 
reference  [4]  for  details).  Table  II  shows  a  comparison  of  the  observed  growth  rate  of  the  sixth-order 
explicit  scheme  with  the  rate  predicted  from  an  eigenvalue  determination. 


Grid 

Q! Numerical 

mm 

usa 

Ffl 

41 

0.1880 

ebsi 

61 

0.1659 

0.1746 

81 

0.1785 

0.1808 

15 


Table  II:  Numerical  vs.  Theoretical  Growth  Rate  for  the  sixth-order  explicit. 


The  agreement  is  very  good,  with  a  slight  discrepancy  in  the  comparison  on  the  61  and  81  grid-point 
cases. 

The  time-divergence  seen  in  the  sixth-order  scheme  is  the  same  as  that  predicted  in  the  counter¬ 
example  presented  in  section  1.  Specifically,  numerical  time  stability  is  not  guaranteed  by  a  dis¬ 
cretization  which  satisfies  a  summation-by-parts  property.  Very  specific  boundary  treatments  must 
be  used  to  guarantee  time  stability. 

SAT  Boundary  Conditions  (Scalar) 

The  SAT  method  for  treating  the  boundary  conditions  guarantees  time  stability  for  the  hyperbolic 
system.  This  method  relies  on  a  spatial  operator  that  satisfies  the  summation-by-parts  energy  norm 
for  the  scalar  case  and  on  very  specific  boundary  treatments  to  ensure  time  stability. 

We  begin  by  showing  that  the  procedure  does  not  destroy  the  formal  accuracy  of  the  spatial 
discretization.  This  result  was  proven  in  section  1  for  the  scalar  case.  Tables  III. a  and  IH.b  show  a 
grid  convergence  study  of  the  SAT  method  on  the  scalar  wave  equation  defined  by  equations  (37), 
(38)  and  (39).  Fourth-order  R-K  time  advancement  is  used  for  all  runs  with  a  time  step  such  that 
no  appreciable  temporal  error  accumulates.  All  calculations  are  run  to  time  T  =  10.  In  all  cases, 
the  calculations  remained  bounded  on  all  grids  (and  CFL’s  less  than  CFLmax)  for  times  as  large  as 
T  =  1000,  which  indicates  time  stability.  This  result  is  consistent  with  the  results  from  eigenvalue 
determinations  in  which  no  RH-P  eigenvalues  were  found. 


r  =  1 

(fourth 

explicit) 

(fourth 

compact) 

(sixth 

explicit) 

Grid 

•og  L2 

Rate 

log  i2 

Rate 

log  L2 

Rate 

mm 

-1.2289 

-1.4005 

-2.5750 

SI 

-2.0878 

4.88 

-2.0479 

3.67 

-3.8300 

7.13 

41 

-2.5784 

3.93 

-2.5096 

3.70 

-4.6500 

6.56 

61 

-3.2211 

3.65 

-3.1689 

3.74 

-5.7880 

6.46 

81 

-3.6806 

3.68 

-3.6464 

3.82 

-6.6056 

6.54 

Table  III. a.  Absolute  error  (logZ/2)  and  convergence  exponent  with  SAT  parameter  r  =  1,  for  the 
fourth  explicit,  fourth  compact  and  sixth  explicit  spatial  discretizations. 
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B 

(fourth 

explicit) 

(fourth 

compact) 

(sixth 

explicit) 

Grid 

log  L2 

Rate 

log  L2 

Rate 

log  L2 

Rate 

wm: 

-1.3472 

-1.8061 

-2.7007 

m 

-2.0866 

-2.4296 

3.54 

-3.8229 

6.37 

41 

-2.5980 

mt 

-2.8773 

3.58 

-4.6666 

6.75 

61 

-3.3107 

MEM 

-3.5243 

3.67 

-5.8518 

6.73 

81 

-3.8145 

4.03 

-3.9978 

3.79 

-6.6485 

6.38 

Table  Ill.b:  Absolute  error  (log  L2)  and  convergence  exponent  with  SAT  parameter  r  =  2,  for  the 
fourth  explicit,  fourth  compact  and  sixth  explicit  spatial  discretizations. 


A  comparison  of  the  SAT  grid  refinement  studies  (table  III. a  and  Ill.b)  with  those  from  the  con¬ 
ventional  boundary  treatment  (table  I),  indicates  that  the  formal  accuracy  of  the  spatial  operator  is 
unaffected  by  the  SAT  treatment.  The  proof  of  stability  given  in  section  1  indicated  that  a  sufficient 
condition  for  stability  of  the  scalar  wave  equation  with  the  SAT  method  is  1  <  r.  The  results  shown 
in  tables  III. a  and  Ill.b  indicate  that  the  magnitude  of  the  error  is  dependent  on  the  value  of  the 
parameter  r.  To  optimize  the  value  of  the  parameter  r  for  these  simulations,  the  e-ror  at  T  =  10 
was  studied  as  a  function  of  r.  An  eigenvalue  code  was  then  used  to  determine  the  maximum  CFL 
of  the  scheme  as  a  function  of  r.  The  results  of  this  study  are  shown  in  Table  IV. 


r 

log  L2 

CFL 

3.0 

-3.8220 

1.17 

2.5 

-3.8221 

1.77 

2.0 

-3.8145 

2.07 

1.75 

-3.8038 

2.07 

1.50 

-3.8833 

2.07 

1.25 

-3.7460 

2.07 

1.00 

-3.6806 

2.07 

0.97 

0.0 

Table  IV:  Absolute  error  (log  L2)  and  CFL  for  various  values  of  the  SAT  parameter  r,  for  the  fourth 
explicit  spatial  operator. 

Note  that  a  fairly  sharp  cutoff  at  the  theoretical  value  of  r  =  1  is  observed  for  the  fourth-order 
explicit  spatial  operator.  (Values  of  r  =  0.93  and  r  =  0.99  were  obtained  for  the  fourth-order 
compact  and  sixth-order  explicit  schemes,  respectively.  In  addition,  precise  agreement  was  obtained 
at  the  r  cutoff  between  the  eigenvalue  determination  and  the  numerical  simulation  of  the  scalar  wave 
equation.)  For  the  fourth-order  explicit  spatial  operator,  the  error  decreased  monotonically  with  r. 
which  suggests  that  the  value  of  r  should  be  as  large  as  possible.  Conversely,  the  maximum  CFL 
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that  is  achievable  with  the  fourth-order  R-K  scheme  decreases  dramatically  at  r  =  2.  A  value  of 
t  =  2  was  determined  to  be  optimal  for  these  studies. 


SAT  Boundary  Conditions  (System) 

The  last  part  of  the  validation  study  is  to  verify  that  the  SAT  boundary  procedure  ensures  stability 
for  the  hyperbolic  system.  Equation  (31)  defines  sufficient  conditions  for  time  stability  ^|^|(2  — 
2\Jl  —  |/?||L|)  <  r  <  [^p](2  +  2^/l  —  |/?||L|)  in  terms  of  r  and  the  boundary  coupling  matrices  L 
and  R.  The  test  case  chosen  is  the  hyperbolic  system 


du  du 

dt  dx 
dv  dv 

dt  dx 


0, 

0,  0  <  X  <  1,  t  >  0 


(41) 


u(0,f)  =  au(0,£),  u(l,£)  = /?it(l,£),  t  >  0 


(42) 


u(x,  0)  =  sin  2nx,  v(x,  0)  =  —  sin  2nx,  0  <  x  <  1,  (43) 

The  exact  solution  for  a  =  0  =  1  is 

u(x,  t)  =  sin  27r(ar  —  £),  v(x,  t)  —  —  sin  27r(:r  +  t),  0  <  x  <  1,  t>0  (44) 

The  case  \cn  R\  =  1  is  neutrally  stable  and  provides  an  extremely  severe  test  of  the  time  stability 
of  a  numerical  method.  No  central  difference  scheme  of  an  order  greater  than  two,  is  time-stable  for 
this  system,  in  spite  of  the  fact  that  the  spatial  operator  is  stable  for  the  scalar  case  (a  =  /?  =  0). 
Examples  include  the  (3-4-3)  compact  and  (3, 3-4-3, 3)  explicit  fourth-order  schemes,  and  the  (52,52- 
6-52,52)  sixth-order  scheme  that  is  shown  in  reference  [4]  to  be  time-stable  for  the  scalar  case.  All 
three  schemes  used  in  the  scalar  analysis  (fourth-order  explicit  and  compact  and  sixth-order  explicit), 
that  satisfy  the  summation-by-parts  property  are  not  time-stable.  In  all  cases,  the  discrete  solution 
of  the  system  defined  by  equations  (41)  through  (44)  diverges  as  time  becomes  large.  Grid  refinement 
shows  Lax  stability  and  an  order  property  for  each  scheme,  but  not  time  stability. 

The  scalar  analysis  demonstrates  a  precise  relationship  between  schemes  that  are  time-stable  and 
the  structure  of  the  eigenvalue  spectrum  that  arises  from  the  discretization  matrix.  Precisely,  if 
RH-P  eigenvalues  exist,  then  numerical  divergence  can  be  expected  from  the  numerical  simulation. 
Unfortunately,  this  statement  is  a  function  of  the  CFL  that  is  used  to  advance  the  solution.  (See 
reference  [4].)  Values  of  the  CFL  can  be  chosen  for  which  no  numerical  divergence  is  experienced  with 
an  R-K  time  advancement  scheme;  for  this  reason  testing  the  numerical  stability  of  various  spatial 
operators  for  the  fully  discrete  system  in  time  is  impractical. 
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The  alternative  is  to  use  the  eigenvalue  structure  of  the  semi-discrete  problem  as  the  test  for 
stability.  If  a  spatial  discretization  operator  has  no  RH-P  eigenvalues,  then  it  is  assumed  to  be  time- 
stable.  A  derivation  of  the  discretization  matrix  operators  for  the  model  hyperbolic  system  [equations 
(41)  and  (42)]  is  presented  in  Appendix  III.  In  addition,  the  structure  of  the  eigenvalues  is  derived. 

For  our  test  system,  we  take  a  =  /3  in  (44)  and  thus  the  sufficient  condition  for  stability  becomes 
(H/].z.°i)  <  r  <  (242^j~a-).  Given  a  value  of  a  and  a  stable  scheme  incorporating  the  SAT 
boundary  treatment  for  the  system,  there  exist  a  range  in  r  for  which  the  time  discretization  is 
stable.  As  in  the  scalar  case,  good  agreement  exists  between  the  theoretical  and  numerical  stability 
limit.  Therefore,  the  agreement  between  the  theoretical  prediction  and  the  numerical  eigenvalue 
determination  was  used  as  a  test  of  the  validity  of  the  theory. 

Table  V  compares  the  stability  limits  of  the  three  high  order  schemes  for  various  values  of  the 
parameter  o;  the  theoretical  limit  is  compared  with  that  predicted  from  the  eigenvalue  determination 
for  the  2x2  system.  The  number  of  grid  points  used  in  each  case  was  101.  A  study  with  61  points 
showed  similar  results.  In  the  study,  is  the  theoretical  value  of  r  based  on  2~2^~t>a  =  r,  and  r ^ 
is  the  value  as  determined  from  the  eigenvalue  determination.  Specifically,  r/v  was  the  smallest  value 
of  r  for  which  the  numerical  eigenvalues  all  had  negative  real  parts.  In  all  cases  the  agreement  was 
very  good,  which  suggests  the  validity  of  the  theory. 


O' 

1.0 

0.99 

0.90 

0.80 

0.50 

Exact 

TT 

2.0 

1.75 

1.39 

1.25 

1.07 

fourth  explicit 

tn 

2.0 

1.75 

mi 

mi 

m 

fourth  compact 

Tn 

2.0 

1.75 

1.39 

1.25 

m 

sixth  explicit 

Tn 

2.0 

1.72 

1.25 

1.01 

1.00 

Table  V:  The  theoretical  and  numerical  stability  limits  of  SAT  boundary  scheme  for  various  values 
of  a. 

In  these  simple  examples,  we  have  demonstrated  that  the  SAT  boundary  procedure  retains  the 
formal  accuracy  of  the  underlying  spatial  operator  and  provides  a  mechanism  to  stabilize  those  spatial 
operators  that  satisfy  a  summation-by-parts  energy  property.  The  resulting  scheme  is  time-stable  for 
both  the  scalar  and  system  case.  The  numerically  predicted  stability  boundaries  for  the  parameter  r 
closely  match  the  theoretical  predictions.  From  a  practical  perspective,  the  numerical  stability  and 
CFL  of  the  fully  discrete  algorithm  are  functions  of  the  value  of  r.  The  choice  t  —  2  seems  to  be 
well  suited  for  both  the  scalar  and  system  cases  and  guarantees  stability  even  for  the  neutrally  stable 
system  case  where  a  =  /?  =  1 . 
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•1.  ( 'ondnsions 


In  this  paper  we  studied  the  stability  and  time  stability  of  the  semi-discrete  hyperbolic  system 
of  partial  differential  equations.  The  spatial  discretizations  considered  were  high  order  (explicit  and 
compact),  and  their  boundary  terms  were  constructed  such  that  the  derivative  matrix  satisfied  a 
summation-by-parts  formula. 

The  following  results  were  obtained: 

1.  A  systematic  way  was  developed  to  obtain  high-order  accurate  derivative  matrices  (includ¬ 
ing  boundary  terms)  having  a  summation-by-parts  property.  The  method  is  illustrated  In- 
finding  explicit  forms  in  the  4th  order  compact  case. 

2.  The  summation- by-parts  property  does  not,  by  itself,  guarantee  the  stability  and  time 
stability  of  the  scheme,  not  even  in  the  scalar  case.  (Refer  to  the  explicit  sixth-order 
example  cited  in  the  text.) 

To  overcome  this  difficulty  we  introduce  the  simultaneous  approximation  term  (SAT)  in 
order  to  account  for  the  effect  of  the  coupling  of  the  physical  boundary  conditions.  The 
SAT  contains  a  free  parameter  r. 

4.  We  give  bounds  on  r  such  that  the  resulting  scheme  for  the  system  (or  scalar)  case,  we 
have  stability  as  well  as  time  stability. 

5.  Numerical  studies  verify  the  theory. 
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APPENDIX  I 


Construction  of  the  Fourth-Order  Compact  Scheme 
We  begin  with  the  semi-discrete  equation  ut  =  A*  u  where  u  —  («i ,  u2, upj)T,  which  results  from  a 
particular  discretization  of  the  equation  ut  =  ux.  The  matrix  A *  is  then  decomposed  as  A*  =  P~ 1  Q. 
The  interior  scheme  used  is  the  fourth-order  compact  scheme  defined  implicitly  as 

1  dui- 1  dui  1  du{+ 1  3  .  ,  /at  i  \ 


1  du,_  i  du{  1  du{+ 1 

4  dx  dx  4  dx 


"(^<+1  Hi— l) 


(AI.  1) 


Note  that  the  interior  scheme  satisfies  the  summation-by-parts  energy  norm  (as  well  as  the  generalized 
norm).  The  matrices  P  and  Q  can  be  written  in  general  form,  with  boundary  closures  of  arbitrary 
size  N  as 

r  „  „  nl  f  9o,o  •  •  •  qo,N  0  ] 


P  = 


PN,N  4 

I  1  I 

4  4 


;  Q  = 


(]N,N  4 

=1  0  2 

4  u  4 


with  the  H  matrix  written  as 


/io,o  •  •  •  fio  ,N 


H  = 


hiw,o  • 


X 

x  y  x 


To  simplify  the  matrix  algebra,  the  following  new  matrices  are  introduced: 


S  = 


0  1  0 

-1  0  1 

0-101 


;  c  = 


'410 

i  1  4  1 

4  0  14  1 


D  = 


y  x  0 

x  y  x 

0  x  y  x 


A  = 


0 

1  0  . 
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Note  that  S,  C,  and  D  are  M  x  M  matrices,  where  M  is  an  arbitrary  number  that  corresponds  to 
the  number  of  interior  points  in  the  discretization.  The  structure  of  the  matrices  is  tri-diagonal  in 
nature.  The  matrix  A  is  N  x  A/,  and  the  only  non-zero  element  is  a/v,i  =  1- 

Thus,  we  can  write  H ,  P,  and  Q  as 


H  x  A  " 

p 

\A' 

Q 

H  = 

x  AT  D 

;  P  = 

.\*T 

C 

;  Q  = 

s 

where  P,  Q ,  and  H  are  the  N  x  N  submatrices  that  involve  the  unknown  quantities  in  the  matrices 
P,  Q,  and  //,  respectively. 

The  spatial  operator  that  involves  P  and  Q  satisfies  the  generalized  summation-by-parts  energy 
norm  if  a  matrix  H  can  be  found  which  simultaneously  symmetrizes  H  P  and  yields  an  H  Q  matrix 
that  is  nearly  skew  symmetric.  By  defining  W  =  H  P  and  V  =  H  Q,  the  matrices  W  and  V  become 


W 


H  P  +  \AAr  x  AC  +  \H  A 
xAT  P  +  \  DAT  DC  +  I  AT  A 


V 


'  HQ  -  A  AT  x  AS  +  \H  A 
mxATQ  -  \  DAT  DS  +  ^  AT  A 


Thus,  the  matrices  W  and  V  are  important  to  the  stability  properties  of  the  spatial  operator. 
Several  notes  about  the  structure  of  W  and  V  should  be  made  at  this  point.  First,  the  matrices 
A  AT  and  AT  A  are  zero  except  for  the  (N,N)  and  (0,0)  elements,  respectively.  Second,  the  matrix 
DC  +  f  AT  A  is  automatically  symmetric  and  it  has  the  same  tri-diagonal  structure  as  the  D  and 
C  matrices.  Third,  the  matrix  D  S  +  ^f  AT  A  is  automatically  skew-symmetric  which  includes  the 
zero  at  the  (0,0)  position.  The  fourth  quadrant  of  W  and  V  automatically  satisfy  the  conditions  on 
the  generalized  summation-by-parts  energy  norm.  The  remaining  conditions  that  W  and  V  must 
satisfy,  written  in  terms  of  the  submatrices  H ,  P,  Q,  C,  D,  S,  and  A,  are 

H  P  =  (H  P)T  (Al.  2) 
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(AI.  3) 


\at  Ht  +  xCT  At  =  -D  AT  +  x  AT  P 
4  4 


HQ  +  (HQ)T  =  '—AAt  +  \60,oI 


(AI.  4) 


-at  ht  +  x  st  at  =  -d  at  -  x  at  q 
4  4 


(AI.  5) 


where  A  <50)o  is  the  non-zero  element  that  occurs  in  the  first  row  and  column  of  the  matrix.  This 
contribution  to  equation  (AI.  4)  allows  for  a  non-zero  value  at  the  (0,0)  element  in  the  matrix  V . 

By  expanding  the  specific  terms  in  equations  (AI.  2  through  AI.  5),  we  ha^e 

^0,n  ■  •  •  hn,n  Pu,0  •  •  •  Pn,n 

0  0  0  0 
AT  HT  =  .  ;  AT  P  = 


ATQ  = 


;  CtAt  = 


•  0  y 


ST  AT 


D  AT  = 


By  comparing  the  matrices  involved  in  equation  (AI.  3),  it  is  apparent  that 

1  y 

-  hk,N  +  x  &k,N  =  xpN,k  +  -  k  =  0,N  (AI.  6) 

Similarly,  equation  (AI.  5)  yields  the  expression 

3  3  y 

^  hk,N  =  —xqNtk  +  —  &k,N\  k  =  0,  N  (AI.  7) 

Eliminating  hk,N  between  equation  (AI.  6)  and  equation  (AI.  7)  yields  the  expression 


<lN,k  =  —3  p/v,fc  +  3  Sk,N,  k  =  0,  N 
24 


(AI.  8) 


These  properties  of  the  matrices  P  and  Q  must  be  satisfied  regardless  of  the  order  properties  of  the 
boundary. 

We  now  derive  the  additional  constraints  that  must  be  satisfied  near  the  boundaries  to  guarantee' 
the  order  properties  of  these  points.  Substitution  of  the  equations  u:  =  jT  and  ^  =  r  jr~'  into 
the  matrices  Q  and  P,  respectively,  yields  the  constraints  that  ensure  the  accuracy  of  the  boundary 
points.  The  general  expression  at  the  boundary  written  in  terms  of  an  arbitrary  accuracy  r  becomes 


N  /Vo 

rY,P^jr~'  +  7  $k,N  (A  +  h,N  (A  +  l)r;  A-  =  0..V  (AI.9) 

J= 0  4  ]=0  4 

Third-order  accuracy  at  the  boundary  points  requires  r  =  0,3  with  A  >  3. 

Thus  far,  we  have  not  specified  the  exact  value  of  the  parameter  A.  We  now  specify  a  precise 
value  for  the  parameter  A  so  that  specific  boundary  conditions  can  be  derived  for  the  fourth-order 
interior  Pade  scheme.  To  retain  the  formal  accuracy  of  the  interior  scheme,  the  boundary  closure 
must  be  accomplished  to  at  least  third-order  accuracy,  and  reqTv.  Tat  A  >  3.  For  .V  =  3. 
equation  (AI.  9)  can  be  written  concisely  in  matrix  notation  as 

‘0*0-'  1*0°  2*0’  3*02]  [0°  01  02  03  ■ 

£  O  +  l"1  1*1°  2*1!  3*12  A  1°  l1  l2  l3 

F  0  *  2"1  \*  2°  2  *  21  3  *  22  ^  2°  21  22  23  + 

.  0  *  3"1  1  *  3°  2  *  31  3  *  32  u  3°  31  32  33  . 

Solving  this  expression  for  the  matrix  Q  results  in  the  expression 

+ 


which  relates  the  matrix  Q  to  the  matrix  P  through  third-order  accuracy  constraints. 

We  will  now  solve  for  the  last  row  of  the  matrices  P  and  Q  and  for  the  last  column  of  the  matrix 
H.  Equation  (AI.  9)  is  written  for  k  =  A,  and  and  (defined  in  equation  (AI.  S))  are  used 
to  yield  the  relationship 


0  0  0  0' 

0  0  0  0 

0  0  0  0 

t_  ^5  rr  ^<23 

24  4  8  12  - 


’0  0  0  O' 

0  0  0  0 
0  0  0  0 
ii  T  10  36 

'-4  4  J 
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r  =  0,3  (AI.  10) 


r  S  PNjJr  1  +  7  +  O"  1 

j— o  4 


—3  ^  Pn,j  jr  +  t  (N  +  l)r  +  3  NT; 

3=0  4 


Setting  N  =  3  and  solving  the  system  for  p3^,  k  =  0,3  yields  p30  =  p3,i  =  0,  p3;i  —  t,  and  p33  =  1. 
Equation  (AI.  8)  can  be  used  to  show  that  <73to  =  ?3,i  =  0,</3|2  =  —  f,  and  q33  =  0.  Similarly, 
equation  (AI.  6)  yields  the  values  of  h^3  as  /i0, 3  =  h\t3  =  0,  h2,3  =  x ,  and  /i3,3  =  y.  Thus,  the  last 
row  of  P  and  Q  are  the  same  as  the  interior  scheme.  In  addition,  the  specific  form  of  the  matrix  H 
must  be 


ho,o 

^0,1 

ho;2 

0  ‘ 

hi,o 

^1,1 

^1,2 

0 

h  2,0 

^2,1 

^2,2 

X 

„  ^3,0 

^3,1 

^3,2 

y  . 

Thus,  accuracy  constraints  on  the  last  row  of  the  matrices  P  and  Q ,  combined  with  the  structure 
requirements  imposed  by  equations  (AI.  3)  and  (AI.  5),  allow  for  the  direct  solution  of  the  last  rows 

A  A  A  A  A  A 

of  P  and  Q,  and  the  last  column  of  H.  Multiplying  the  expression  relating  P  to  Q  by  the  matrix  H, 
and  using  the  substitutions  H  Q  =  V  and  H  P  —  W  yields  the  expression  for  V  of  the  form 


V  = 


3 

-1 

2 

-1 

3 
2 


‘0  0  0  0‘ 

0  0  0  0 

7  x  —5  x  17  x  —23  x 

L  24  4  8  12  J 


Solving  for  W  and  V  such  that  equation  (AI.  2)  (where  W  =  WT)  and  equation  (AI.  4)  (V  +  VT  = 
A  AT  +  A  60,o  /)  are  satisfied  to  obtain 


[  -9  a 

16 

1536  7+1536/3- 899  a 

768 

768  7+768/3- 703  a 
192 

15367+1536/3-1481  a 
768 

15367+1536/3-899  a 

0 

1536  7+1536 /3-1277a 

7687+768/3-733  a 

768 

256 

192 

768 -y+768 /3— 703  a 

1536t(+1536  0—1277  a 

0 

15367+1536/3-947a 

192 

256 

768 

1 536-/+ 1536/0 -1481  a 

768-^+768/3—733  a 

15367+1536  0-947  a 

-3  a 

L  768 

192 

768 

32 

and 
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w  = 


137  ft  — 192  -y 
192 

—512  "y— 128  /3+525  a 
128 

_ 145  ft— 144  -y 

48 


-512 -y- 128/3+525  ft 
128 


1 45  ft—  1 44  *y 
48 


— 816~y— 384  /3+737  a  -3072  -y- 1344  /3+3001  a 

48  192 

-3072^-1344/3+3001  ft  -3264  -y-1536  /3+301 1  « 


P 


192 

557  ft  — 576  -y 
192 


192 

-1536  ->-384/3+1543  0 
384 


P 

557  ft  — 576  -y 
192 

-1  536  ft  — 384  /3+1543 ft 
384 


with  x  =  -jp  and  y  =  a.  Three  arbitrary  parameters  remain  after  all  accuracy,  symmetry  and 
skew-symmetry  conditions  are  satisfied. 

The  final  step  in  the  discretization  is  to  find  a  specific  form  of  the  matrix  P  that  will  lead  to  a 
simple  algorithm.  Because  the  matrix  P  is  tri-diagonal  in  the  interior,  the  boundary  closure  should 
retain  the  tri-diagonal  structure.  After  P  is  specified,  we  can  solve  for  the  matrix  H  from  H  =  V  P~x 
if  the  inverse  of  P  exists,  and  the  last  column  of  H  is  [0, 0,  y,  x]T.  The  matrix  Q  follows  immediately 
from  Q  —  P  V~x  W .  The  last  test  is  to  ensure  that  both  W  and  that  the  full  matrix  W  are  positive 
definite. 

Many  matrices  P  have  been  found  that  satisfy  all  of  the  criteria  given  in  the  generalized  summation- 
by-parts  energy  norm  analysis.  From  a  numerical  perspective,  all  behaved  similarly.  The  results 
presented  here  are  those  that  were  the  simplest  to  code.  Choosing  a  specific  matrix  P  of  the  form 


P  = 


yields  a  matrix  Q  of  the  form 


-289 

279 

75 

—  7 

234 

286 

286 

2574 

-8635 

6987 

1851 

-203 

3376 

3376 

3376 

3376 

-15043 

-4089 

147 

29353 

18972 

2108 

124 

18972 

0 

0 

-3 

A 

0 

r  m  1  0  0 

429  1  v  v 

1  3563  r, 

1  1688  8  U 

n  43  1893  139 

U  17  1054  186 

0  0  \  1 

4 


The  resulting  matrix  H  is  therefore 
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70282007653 

7658388480 

-9426299 

2268480 

-192913 

1067520 
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-55530689643 

2552796160 

8051589 

756160 

1 49823 
355840 
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63842626133 

2552796160 

-9153739 
7561 60 

-4433 

355840 

-1 

8 

-71498870443 

7658388480 

10110149 

2268480 

102703 

1067520 

1 

^■.Froin  a  practical  point  of  view,  the  inconvenient  form  of  the  H  matrix  is  not  of  great  concern 
since  the  matrix  H  is  only  inverted  once  and  one  column  is  stored  for  use. 

The  matrices  P,  Q ,  and  H  can  be  used  to  establish  both  the  symmetry  of  the  matrix  V  and  the 
near  skew-symmetry  of  the  matrix  W .  The  first  six  rows  and  columns  of  the  V  matrix  are 


16513 

-261 

2993 

-6223 

0 

46080 

5120 

15360 

46080 

-261 

9153 

-2943 

1611 

0 

5120 

5120 

5120 

5120 

2993 

-2943 

7473 

-2063 

-1 

15360 

5120 
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-1 

32 
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16 

0  0  0  —  -  —  I 
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The  first  six  rows  and  columns  of  the  W  matrix  are 
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0 
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32 

-3 

4 

0 

3 
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0 
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0 

3 

32 

-3 

4 
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As  shown,  the  matrix  W  is  nearly  skew  symmetric,  and  the  matrix  V  is  symmetric.  For  the  matrix  IF 
is  positive  definite,  it  is  necessary  to  show  that  every  submatrix  is  positive  definite.  The  inner  scheme 
is  diagonally  dominant  and  contributes  to  the  definiteness  of  the  complete  matrix  W.  However,  the 
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boundary  elements  are  not  diagonally  dominant,  and  suppress  the  positive-definiteness.  The  l  x  1 
boundary  matrix  W'  —  VV'  -f-  |  A  A1  has  the  following  characteristic  polynomial 

754974720  A4  -  3507814400  A3  +  5299068928  A2  -  3079323424  A  +  530779791  =  0  (Al.  11) 

The  symmetry  of  the  IT'  matrix  and  the  alternating  signs  of  each  term  in  the  characteristic  polynomial 
guarantee  that  the  matrix  is  positive  definite.  The  characteristic  polynomial  of  every  submatrix  (up 
to  ten  points,  which  includes  four  boundary  and  six  interior  points)  of  the  matrix  It  results  in  a 
positive  definite  matrix.  No  proof  that  the  complete  discretization  is  positive  definite  for  an  arbitrary 
number  of  interior  points  has  been  found. 

The  accuracy  of  the  new  scheme  is  third  order  at  the  boundaries  anil  fourth  order  in  the  interior. 
To  show  this,  the  Taylor  expansion  for  long  wavelength  modes  is  made  using  the  stencil  at  each  of 
the  first  four  points.  The  results  are 


u  +  +  "' 
*... 

^  2005 5 


*  4 


i - r  + 

18(T 


(Al.  12) 


At  high  resolution,  the  boundary  points  behave  with  third-order  truncation  error;  the  interior  behaves 
with  fourth-order  error.  Therefore,  the  resulting  scheme  is  formally  fourth-order  accurate. 
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APPENDIX  II 


■Sixth  Order  Explicit  Scheme 

Here,  we  derive  an  explicit  scheme  that  is  formally  sixth-order  accurate.  Unlike  the  fourth-order 
compact  case  presented  earlier,  the  matrix  H  can  be  the  identity  matrix.  To  constrain  the  matrix  P 
to  be  symmetric  and  the  matrix  Q  be  nearly  skew  symmetric,  six  alternative  formulas  are  required 
at  the  boundaries,  each  of  which  is  closed  to  fifth-order  accuracy  to  retain  the  formal  accuracy.  The 
corner  7x7  submatrices  of  the  global  matrices  P  and  Q  can  be  written  as 
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The  characteristic  polynomial  of  the  matrix  Pe,  is 


1 0399739562845798400000000  A6 
+  1003578630643249838161920000  A4 


2485 1 26099 1 6244983808000000  A5 
1 63903822337723736805 1 7 1 2000  A3 


+  1248376737213799711434406800  A2  -  4 1 22353650428 1 6633559 1 97440  A 

+  37455444120716264727507839  =  0  (All.  1) 

The  symmetry  of  the  matrix  P$  and  the  alternating  signs  of  the  terms  in  the  polynomial  are  sufficient 
for  positive  definiteness  of  both  the  matrix  P$  and  the  global  matrix  P. 

The  truncation  error  at  the  boundary  points  is 


*£  + 
it,  - 
it  - 
it  + 
it  + 
it  ~ 


644829999745 1 547397244467^ 
224732664724297588365047034 
55 1 7845934 1 997062554732 1 £® 

1 1 2366332362 1 48794 1 825235 1 70 
90378 1 1 40428 1 60989627296 1 9£® 
2247326647242975883650470340 
62520732887440126777806839^® 
2247326647242975883650470340 
2 1 52 1 02 1 0826949659 17733 1 
1 1 2366332362 1 48794 1 825235 1 70 
7101 580254 19711 6302053905^® 
224732664724297588365047034 


(All.  2) 


which  indicates  fifth-order  accuracy  at  the  six  boundary  points. 
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APPENDIX  III 


Eigenvalues  of  the  Discrete  System 

The  eigenvalues  of  the  semi-discrete  system  are  used  in  the  results  section  to  compare  the  theoretical 
and  the  numerical  stability  boundaries.  The  model  equation  is  the  hyperbolic  system  used  in  the 
main  text  and  defined  by  equations  (41),  and  (42).  For  convenience,  we  define  the  (.V  +  1 )  x  (.V  -f  1) 
matrix  A  =  P_1  Q.  The  matrix  A  contains  all  the  information  from  the  spatial  discretization 

operator  The  semi-discrete  form  of  equation  (41)  becomes 


du 

Tt+Au 


dv 

dt 


—  A  v 


0. 

0,  t  >  0 


( A 1 1 1 .  1) 


with  the  boundary  conditions  defined  by  equation  (42).  In  matrix  notation,  the  discrete  system  takes 
the  form 
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Note  that  JJ  =  /,  so  that  J  —  J-1.  The  vector  u  is  the  concatenated  vector  of  discrete  values 
from  the  scalar  vectors  u  and  v  with  the  elements  Uo  and  v\<  removed.  These  elements  are  removed 
because  the  physical  boundary  condition  relates  them  to  known  elements  in  the  vector  u,  so  that  and 
need  not  need  to  solve  for  them.  The  matrix  A *  is  the  N  x  N  submatrix  of  A  which  is  obtained  by 
eliminating  the  zeroth  row  and  zeroth  column.  Note  that  this  was  the  matrix  that  was  analyzed  in  the 
scalar  analysis  to  determine  time  stability  of  the  spatial  operator.  The  matrix  B  is  zero  everywhere 
except  the  first  column,  where  the  zeroth  column  of  the  original  A  matrix  is  written.  This  column  is 
precisely  the  coupling  between  the  u  and  v  vector,  which  occurs  at  the  boundary. 

It  is  instructive  to  relate  the  system  eigenvalues  to  those  obtained  in  the  scalar  analysis  [(/l*  — 
X  I)u  =  0].  By  defining  the  matrix  H~l  and  H  as 


'  y/Bl 

\/ct  J 

\/otl  . 

—  y/ttl 

H~ 1  =  -yd— 

V2v^  [  -s/B  I 

y/a  J 

h  =  -y—  .  . 

V2V^  [  s/BJ  . 

with  H  1  H  =  H  H  1  =  /,  we  note  that  the  system  matrix  can  be  made  block  diagonal  with  the 
similarity  transform  H 


1 

V0I 

.  y/ctj  ' 

'  aB  ' 

'  y/Sl  . 

->/  1 

i 

\J 2  l\/o0 

-VBi  ■ 

■  \/aJ 

(3J~lBJ  .  J-MtJ 

.SJ  ■ 

yfflj  V2V^ 

'  A f  +  yjcrf  BJ  .  0 

0  .  -  s/7^BJ 

For  scalar  time-stable  spatial  schemes,  the  eigenvalues  of  the  matrix  A t  are  bounded  to  the  left  half¬ 
plane.  Note  that  for  a  =  0  (or  {3  =  0)  the  contribution  from  the  boundary  coupling  matrix  B  is 
identically  zero,  and  the  eigenvalues  of  the  resulting  system  are  simply  the  scalar  eigenvalues  with 
a  multiplicity  of  two.  For  non-zero  values  of  the  parameters  a  and  /?,  the  eigenvalues  of  the  total 
matrix  are  different  from  those  of  the  original  matrix  A h  Also  note  that  two  distinct  eigenvalue 
scenarios  exist  for  the  boundary  parameters  «  and  (3,  depending  on  whether  their  signs  are  equal  or 
opposite. 
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